Introduction

Here we will recreate the ZooMSS model (version 2) in Heneghan et al. (2020 in preparation) using mizer.

We begin with some setup of required packages.

#get required packages
library(devtools)
#most up to date master branch of mizer
#install_github("sizespectrum/mizer")
#install_github("astaaudzi/mizer-rewiring", ref = "temp-model-comp")
#documentation here:
#https://sizespectrum.org/mizer/dev/index.html
library(mizer)
require(tidyverse)

#remotes::install_github("sizespectrum/mizerExperimental")
library(mizerExperimental)

Data

Next, let’s load in and check the data used in the ZooMSS model. TODO: add in some biomass data

#read in group data
groups <-read.csv("data/TestGroups_mizer.csv")
groups$w_min <- 10^groups$w_min
groups$w_inf <- 10^groups$w_inf
groups$w_mat <- 10^groups$w_mat
#groups<-merge(groups,dat,by.x="species",by.y="group",all=T)

# have a look at plot
# plot(groups$w_inf,groups$biomass.tperkm2,xlab="log Maximum Weight [g]", ylab=" log Total Biomass", log="xy",col="blue",pch=16)
# text(groups$w_inf,groups$biomass.tperkm2,labels=groups$species,cex=0.5)

# could plot the paramter allometries here to explore

Set-up mizer model

Next let’s read in the parameters from ZooMSS.

#read groups again for southern ocean model, this time subsetting key groups
groups <-read.csv("data/TestGroups_mizer.csv")
groups$w_min <- 10^groups$w_min #convert from log10 values
groups$w_inf <- 10^groups$w_inf
groups$w_mat <- 10^groups$w_mat
groups$h <- 10000000000 # should be Inf, but breaks. Massive value still works out to effectively unlimited feeding as allowed in ZooMSS

#groups <- readRDS("data/groups.RDS")[-1,]
#check fails these tests:
#groups$w_mat25 >= groups$w_mat #note we don't have w_mat25. Not sure if it's necessary/desirable. Probably not.
groups$w_mat <= groups$w_min
##  [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
groups$w_inf <= groups$w_mat
##  [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
# read interaction matrix
# get the interaction matrix - actually I think we can leave this out. Default is all 1s, which is the same as in ZooMSS. Included for completeness; it may be useful in future to keep this in.
theta <- readRDS("data/zoomss_inter.rds")[,-1]
#[-1,-1]

We will pass these parameters to mizer to set up a new multispecies model.

TODO: adjust parameters here.

ID <- 1 #index of environmental data to choose
enviro_row <- readRDS("data/envirofull_20200317.RDS")[ID,]

#set up the fixed phyoplankton spectrum
phyto_fixed <- function(params, n, n_pp, n_other, rates, dt = 0.01, phyto_int = enviro_row$phyto_int, phyto_slope = enviro_row$phyto_slope,...) {
    n_pp <- 10^(phyto_int)*(params@w_full^(phyto_slope)) #returns the fixed spectrum at every time step
    n_pp[params@w_full>params@resource_params$w_pp_cutoff] <- 0
    return(n_pp)
}

mf.params <- newMultispeciesParams(species_params=groups,
                                   interaction=NULL, #NULL sets all to 1
                                   no_w = 178, #number of zoo+fish size classes;
                                   kappa = 1e15, #set huge to recreate unbounded feeding on phyto in ZooMSS
                                   min_w_pp = 10^(-14.4), #minimum phyto size. Note use -14.4, not -14.5, otherwise it makes an extra size class
                                   w_pp_cutoff = 10^(enviro_row$phyto_max), #maximum phyto size
                                   r_pp = 20, #Coefficient of the intrinsic resource birth rate, set large for now
                                   n = 0.7, #The allometric growth exponent used in ZooMSS
                                   z0pre = 1, #external mortality (senescence)
                                   z0exp = 0.3,
                                   resource_dynamics = "phyto_fixed",
                                   RDD = constantRDD(species_params = groups) #first go at this
                                   #pred_kernel = ... #probably easiest to just import this/pre-calculate it, once dimensions are worked out
)
## Note: No ks column so calculating from critical feeding level.
## Note: Using z0 = z0pre * w_inf ^ z0exp for missing z0 values.
#checks
length(which(mf.params@initial_n_pp>0)) == length(seq(-14.5,enviro_row$phyto_max, by = 0.1))
## [1] TRUE

Now do some fiddling to make the new MizerParams object match the ZooMSS parameters. Note: this chunk is adapted from fZooMSS_setup.R, found at https://github.com/MathMarEcol/ZooMSS/.

setZooMizerConstants <- function(params, Groups, input_params){
#### CALCULATES CONSTANT BITS OF THE MODEL FUNCTIONS FOR EACH GROUP
  SearchVol <- getSearchVolume(params)
  M_sb <- getExtMort(mf.params)
  ZSpre <- 1 # senescence mortality prefactor
  ZSexp <- 0.3 # senescence mortality exponent
  
  pred_kernel <- getPredKernel(params)
  prey_weight_matrix <- matrix(params@w_full, nrow = length(params@w), ncol = length(params@w_full), byrow = TRUE)
  pred_weight_matrix <- matrix(params@w, nrow = length(params@w), ncol = length(params@w_full))
  
  for(i in 1:nrow(params@species_params)){
    ## Senescence mortality
    if(params@species_params$Type[i] == "Zooplankton"){
      M_sb[i,] <- ZSpre*(params@w/(10^(params@species_params$w_mat[i])))^ZSexp
      M_sb[i, 10^(params@species_params$w_max[i]) < params@w] <- 0
      M_sb[i, 10^(params@species_params$w_mat[i]) > params@w] <- 0
    }
    
    if(params@species_params$Type[i] == "Fish"){
      M_sb[i,] <- 0.1*ZSpre*(params@w/(10^(params@species_params$w_mat[i])))^ZSexp
      M_sb[i, 10^(params@species_params$w_max[i]) < params@w] <- 0
      M_sb[i, 10^(params@species_params$w_mat[i]) > params@w] <- 0
    }
    
    ### Search volume
    SearchVol[i,] <- (params@species_params$gamma[i])*(params@w^(params@species_params$gamma[i]))
    SearchVol[i, 10^(params@species_params$w_max[i]) < params@w] <- 0
    SearchVol[i, 10^(params@species_params$w_min[i]) > params@w] <- 0
    
    ### Predation Kernels
    if(is.na(params@species_params$PPMRscale[i]) == FALSE){ # If group has an m-value (zooplankton)
      # Calculate PPMR for zooplankton, which changes according to body-size (Wirtz, 2012)
      D.z <- 2*(3*params@w*1e12/(4*pi))^(1/3) # convert body mass g to ESD (um)
      betas <- (exp(0.02*log(D.z)^2 - params@species_params$PPMRscale[i] + 1.832))^3 # Wirtz's equation
      beta_mat <- matrix(betas, nrow = length(params@w), ncol = length(mf.params@w_full))
      
      # Calculate feeding kernels
      pred_kernel[i, , ] <- exp(-0.5*(log((beta_mat*prey_weight_matrix)/
                                     pred_weight_matrix)/params@species_params$FeedWidth[i])^2)/
        sqrt(2*pi*params@species_params$FeedWidth[i]^2)
      # The feeding kernel of filter feeders is not expected to change much with increasing size so we fix it here
      
      # if (param$fixed_filterPPMR == TRUE){
      if(i == 3){
        pred_kernel[i, , ] <- matrix(pred_kernel[i,44,], nrow = length(params@w), ncol = length(mf.params@w_full), byrow = TRUE)
      }
      if(i == 8){
        pred_kernel[i, , ] <- matrix(pred_kernel[i,61,], nrow = length(params@w), ncol = length(mf.params@w_full), byrow = TRUE)
      }
      # }
      
    } else { # If group does not have an m-value (fish)
      beta_mat <- matrix(params@species_params$PPMR[i], nrow = length(params@w), ncol = length(params@w_full))
      
      # Calculate feeding kernels
      pred_kernel[i, , ] <- exp(-0.5*(log((beta_mat*prey_weight_matrix)/
                                             pred_weight_matrix)/params@species_params$FeedWidth[i])^2)/
        sqrt(2*pi*params@species_params$FeedWidth[i]^2)
    }

  }
  
  #temperature effect 
  temp_eff <-  matrix(2.^((input_params$sst - 30)/10), nrow = length(params@species_params$species), ncol = length(params@w))
  M_sb <- temp_eff * M_sb # Incorporate temp effect on senscence mortality
  
params <- setExtMort(params, z0 = M_sb)
params <- setSearchVolume(params, SearchVol)
params <- setPredKernel(params, pred_kernel)
}

setZooMizerConstants(mf.params, groups, enviro_row)

Try running it:

sim <- project(mf.params, t_max=5,dt = 0.01)
plot(sim)

plotlyBiomass(sim)
#plotlyGrowthCurves(sim,species="macrozooplankton")
plotlyFeedingLevel(sim)
# feeding level satiation for some groups, except for the seabirds
# macrozooplankton - they are not growing enough,why?
#tuneParams(mf.params)
plotlyGrowthCurves(sim,percentage = T)
plotlySpectra(sim)

—–Old stuff only below this line—–